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We perform intensive numerical simulations of the three-dimensional site-diluted Ising antifer- 
romagnet in a magnetic field at high values of the external applied field. Even if data for small 
lattice sizes are compatible with second-order criticality, the critical behavior of the system shows 
a crossover from second-order to first-order behavior for large system sizes, where signals of latent 
heat appear. We propose "apparent" critical exponents for the dependence of some observables with 
the lattice size for a generic (disordered) first-order phase transition. 
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I. INTRODUCTION 

The study of systems with random fields is of 
paramount importance in the arena of the disordered sys- 
tems. A paradigm in this field is the random field Ising 
modelii^ (RFIM). In spite of much effort devoted to the 
investigation of the RFIM;^ several important questions 
remain open. Some of these questions refer to the nature 
of the (replica symmetric?) low temperature phase, to 
universality issues (binary versus Gaussian external mag- 
netic field^i^), and to the order of the phase transitions. 
Here we study the diluted antiferromagnetic Ising model 
in an external magnetic field (DAFF) that is believed 
to belong to the same universality class of the RFIM 
(the DAFF is expected to behave as a Gaussian RFIM 
because of the short ranged correlations in the superex- 
change coupling) i^i^ 

As a matter of fact, DAFF systems are the most widely 
investigated experimental realization of the RFIM. One 
of the best examples of a diluted Ising antiferromagnet 
is Fe:EZni_2.F2. Its large crystal field anisotropy persists 
even when the Fc ions arc diluted (x < 1), thus providing 
a good (antiferromagnetic) Ising system for all ranges of 
the magnetic concentration. Other systems behaving as 
Ising (anti)ferromagnets are Fe2:Mg;^_^Cl2, CoZni_2:F2 
and Mn2;Zni_2;F2 The experimental results on the or- 
der of the phase transition are somewhat inconclusive. 
On the one hand, these materials show a large critical 
slowing down around the critical temperature, as well as 
a symmetric logarithmic divergence of the specific heat. 
On the other hand, the order parameter (the staggered 
magnetization) behaves with a critical exponent (3 near 
zero, possibly marking the onset of a first-order phase 
transition.— Note that (3 should be exactly zero if the or- 
der parameter is discontinuous as a function of tempera- 
ture or the magnetic field. 



The numerical investigations of the DAFF at T > 
are scarce. It was investigated a long time ago by Ogiel- 
ski and Huscj^ They considered several values of the pair 
temperature-magnetic field on lattice sizes up to i = 32 
but far from the critical region^ They investigated the 
thermodynamics as well as the (equilibrium) dynamical 
critical behavior. Their thermodynamic results pointed 
to a second-order phase transition. However, they found 
activated dynamics (which could be interpreted as a sig- 
nal of a first-order phase transition) rather than standard 
critical slowing down (for numerical studies of the DAFF 
at T = 0, see Refs. and [l]). 

Numerical and analytical studies rather focused on the 
RFIM, which is expected to display the DAFF critical be- 
havior.i^ Even if the RFIM is more amenable than the 
DAFF to analytical investigations, the situation is still 
confusing. Indeed, mean-field theory predicts a second- 
order phase transition for low magnetic field. If the prob- 
ability distribution function of the random fields does not 
have a minimum at zero field, the transition is expected 
to remain of the second-order all the way down to zero 
temperature. However, if the probability distribution for 
the random field does show a minimum at zero field, a 
tricritical point and a first-order transition line at suffi- 
ciently high field values are predictedJ^ 

The numerical investigation of the RFIM has neither 
confirmed nor refuted this counterintuitive mean-field 
result. Rieger and Young studied binary distributed 
quenched magnetic fieldsjii where mean-field predicts a 
tricritical point. After extrapolation to the thermody- 
namic limit, they interpreted their results as indicative 
of a first-order transition for all external field strengths 
in the thermodynamic limit (the tricritical point did not 
show up). Rieger studied the case with Gaussian fields,— 
where only second-order behavior should be found ac- 
cording to mean- field expectations. Actually, his results 
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were consistent with a second-order phase transition with 
vanishing (!) order parameter exponent. The simulation 
of Hernandez and Diep*^ of the binary RFIM supports 
the existence of a tricritical point at finite temperature 
and magnetic field. Also the study by Machta et al^ 
of the Gaussian RFIM showed evidences of finite jumps 
in the magnetization at disorder dependent transition 
points. 

A completely different numerical strategy is suggested 
by the expectations of a T = renormalization-group 
fixed point^. Since the ground state for a RFIM on a 
sample of linear size L can be found in polynomial (in L) 
time, T = physics can be directly addressed by study- 
ing the properties of the ground state for a large number 
of samples. Hartmann and Youngi^ studied lattices up 
to L = 96 for the gaussian RFIM. They concluded that 
their data supported a second-order phase transition sce- 
nario for the Gaussian RFIM. The same model was in- 
vestigated by Middleton and Fisher— in lattices up to 
i = 256. Their data suggested as well a continous phase 
transition with a very small order parameter exponent, 
/3 = 0.017(5). Using the same technique Hartmann and 
Nowak^ found non-universal behavior in the binary and 
Gaussian RFIM, not excluding that the former could un- 
dergo a first-order phase transition. They also studied 
T = critical properties of the DAFF for system sizes 
up to L = 120, and found critical exponents /? = 0.02(1), 
1^^:1.14(10) and 7 = 3.4(4), compatible with their results 
for the Gaussian RFIM. 

The aim of this work is to revisit the Ogielski-Huse in- 
vestigation, that was carried out for T > 0, with modern 
computers, algorithms^ and finite-size scaling methods 
(the quotient methodfi^ that uses the finite volume cor- 
relation length2S to characterize the phase transition). 
The significant CPU investment allowed us to simulate 
large lattices (L = 24) and a large number of disorder 
realizations. We plan, in the future, to perform a more 
complete investigation of the critical surface of this model 
(this would require to vary three variables: temperature, 
dilution and magnetic field). 

Our main finding is that the DAFF probably under- 
goes a very weak first-order phase transition. This seems 
a natural explanation for the finding of the activated dy- 
namics (at equilibrium) in Ref. [7|. The discontinuity 
in the magnetization density is sizable, yet that of the 
internal energy is very small. Nevertheless, even if we 
perform a standard second-order analysis, the critical ex- 
ponent for the staggered magnetization turns out to be 
ridiculously small. 

The outline of the rest of this paper is the following: in 
the next section we describe the model (Sect. Ill A[) . ob- 
scrvables (Sect. IiTB|) . as well as the theoretical expecta- 
tions for the finite-size scaling behavior in a second-order 
phase transition (Sect. Ill C[) and for a first-order one 
(Sect. IIID[) . Details about our simulations are given in 
Sect. mil Our numerical results are presented in Scct. lIVI 
We first perform a second-order analysis (Sect. IIV A[) . 
then consider the possibility of a weak first-order tran- 



sition (Sect. IIV Bp . After summarizing our results in 
Sect. |Vl we discuss in the Appendix that the boundSi 
ly > 2/d oi Chayes et al. holds as well for first-order 
phase transitions in the presence of disorder. In addi- 
tion, we have found upper bounds for the divergence of 
the susceptibility and specific heat with the lattice size. 

II. MODEL 

A. Model, phase diagram, symmetries 

The model is defined in terms of Ising spin variables 
Si = ±1, i = 1, . . . , = L"^ placed on the nodes of a cubic 
lattice of size L with periodic boundary condition. The 
spins interact through the following lattice Hamiltonian 

7i = ^ e^SiCjSj - H^^eiSi , (1) 

<i,j> i 

where the first sum runs over pairs of nearest-neighbor 
sites, while H is the external uniform magnetic field. The 
Ci are quenched dilution variables, taking values and 1 
(with probability 1 — p and p respectively) in empty and 
occupied sites. In this work we have fixed this probability 
to p = 0.7. In this way we arc guaranteed to stay away 
from both from the pure case (p = 1) and the percolation 
threshold {pc w 0.31)^ 

It is understood that for every choice of the {ei}Y=i, 
called hereafter a sample or a disorder realization, we are 
to perform the Boltzmann (thermal) average. The mean 
over the disorder is only taken afterwards. 

For low magnetic field, at low temperatures, model 
P]) stays in an antiferromagnetic state that we will call 
the ordered phase. The staggered magnetization Mg, see 
Eq.(I2]) below, is an order parameter for this phase. Note 
that, for iJ = 0, the Z2 transformation 5*^ — > —Si, yields 
a degenerate antiferromagnetic state. The increase of 
the magnetic field or the temperature T weakens the an- 
tiferromagnetic correlations, and the system eventually 
enters into a paramagnetic state. The paramagnetic and 
the ordered phases are separated by a critical line in the 
{T,H) plane [it would be a critical surface in the {T,H ,p) 
phase diagram]. 

Note that the effect of disorder (quenched dilution), 
combined with the applied field H va a. finite DAFF sys- 
tem, breaks the Z2 symmetry even in the ordered phase. 
Consider the state of minimal energy at T = and H low 
enough so that the staggered magnetization is maximal 
and Ms = p. Now let us change the sign to all spins in one 
hit: if p = 1, the two states arc completely degenerate, 
but the random dilution introduces a subextensive shift 
in the total energy. In fact, the inversion does not change 
the nearest-neighbor energy, but changes the sign of the 
magnetic part. In the pure system the magnetic energy 
of the fully ordered antiferromagnetic state is zero but, 
in the presence of random dilution, the number of spins 
aligned or misaligned with the field iJ is a random vari- 
able. So, the total magnetic contribution to energy is of 
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1 /2 

order [p{l — p)N] ' . It follows that the Ms = —p state 

1 /2 

has an energy shift of order 2 [p{l — p)N] with respect 
to the Ms = p one and its Boltzmann weight is depressed 
(an analogous effect of degeneration removal is present in 
the RFIM). A "quasisymmetric" state may exist if there 
is a configuration of spins in which almost all the spins 
are reversed with respect to the Ms = p state such that 
the sum of the energies of unsatisfied bonds cancels the 
magnetic excess. If this is the case, then the two states 
are degenerate but the probability distribution of the or- 
der parameter results peaked around asymmetric values. 
The magnetic energy excess is a subextensive effect and 
is expected to be suppressed as L increases. Nevertheless, 
the probability of transitions between states of opposite 
spontaneous staggered magnetizations decreases for large 
system, which is a major problem for simulations. 



B. Observables 



In the following, ((• • •)) denotes thermal averages (in- 
cluding averages of real replicas) and (• • •) indicates a 
sample average (average on the disorder). Measures fo- 
cused on several observables: the order parameter, i.e., 
the average value of the staggered magnetization 



(2) 



(j^ is the fi-th lattice coordinate of site j) whose values 
are limited in the interval —p < Mg < p (in average for 
large lattices); the average energy densities are 
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(n) = {e) = {Ek) + h{Em), 



(3) 



^{Hk) = {Ek) = ^{J2 ^^S,e,S,), (4) 



<ij> 



(5) 



with Hk and Hm respectively the kinetic and magnetic 
contributions to the Hamiltonian. The definition of Em 
coincides with the definition of the usual magnetization 
density. 

We also computed the average values of the squares 
and fourth powers of the above quantities, and some cu- 
mulants and susceptibilities: given an observable yA — 
A we compute the Binder cumulant: 
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and the connected and disconnected susceptibilities 



= V{A^-{AY 
xL = V{Af. 



(6) 



(7) 
(8) 



These are the ordinary susceptibilities in case A = M^, 
while xV' is proportional to the specific heat Cy . 

The lack of Z2 symmetry, explained in Sec. Ill Al makes 
mandatory the use of connected correlation functions in 
finite lattice sizes, especially in the case of the order pa- 
rameter Mg. Yet, the connected staggered susceptibil- 
ity Xc^" = V{MJ^) — {Ms)^ does not show a peak in the 
(T, H) ranges we considered, so we also study the behav- 
ior of the connected and disconnected staggered suscep- 
tibilities defined with the absolute value of the staggered 
magnetization: 



Xdi 



Xc = V{M^ {\Ms\Y 
V{\Ms\f. 



(9) 
(10) 



In the following, when no observable subscript is specified 
in the susceptibility symbol, we will be referring to Eqs. 
di) and PI). 

It will turn out useful to define a correlation length on 
a finite lattice by the following analogy with a (lattice) 
Gaussian model)^ 



G{ki)-G{k2) 
klG{k2) - kjG{ki) ' 



(11) 



with fc2 ^ 4X;):=iSin2 (fcp/2) on a discrete lattice and 
G{k) the momentum-dependent propagator 



G(fc) = V{F{k)F{~k))-{F{k)){F{-k)), (12) 



1 V- ,E^.i(fc.+-b. 
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(G(0) = xf% F(0)=M,), 



(13) 



and F{k) the staggered Fourier transform of the spin 
field. Also in this case we use the connected part for 
G{k). Choosing ki = (0,0,0) and ^2 = (27r/i)fc^ as one 
of the minimum wave vectors (fc^, /i = 1, . . . , c? are the d 
versors in the reciprocal space), we have 
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4sin^(Vi) \ [j:''^^,G[i2n/L)k,)]/d 



(14) 



Equation (fH)) is a good estimate of the correlation length 
only in the disordered phase but is useful to identify the 
critical region where ^/ L should be a nontrivial universal 
value. 

Finally, with mass storage not being a problem on 
modern equipment, it is easy to compute derivatives with 
respect to inverse temperature P = I /T and applied field 
H, through connected correlations. In particular, the 
specific heat is 



Cy 



1 d{n) 

V dT 



(15) 
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Finite-size scaling in second-order phase 
transitions 



D. Finite-size scaling for first-order phase 
transitions 



We made use of finite size scaling,— both studying 
the behavior of peaks of susceptibilities and applying the 
quotient method (QM)i^ to extract values for critical ex- 
ponents. Let us briefly recall both. 

Consider an observable A, that in the infinite volume 
limit behaves as (T — Tc)^" — t^"^ near the critical re- 
gion (t is the reduced temperature). Then, disregarding 
corrcction-to-scaling terms, wc expect the following tem- 
perature dependency on a finite lattice of size L 



A{L,t)=L-'/-fA{tL^'n: 



(16) 



with V the correlation length exponent, ^ oc t^"^, and 
fA{s) a smooth universal scaling function showing a peak 
at some value s™ = tm{L)V'/^ . It follows that Tm{L) — 
oc L~^/^ . In addition, the scaling of the peak- height 
gives the value of the critical exponent a. 

The QM is based on the same scaling ansatz: 



(17) 



We compare data in two lattices Li and L2 at the 
(unique) reduced temperature t* where the correlation- 
length in units of the lattice size coincides, ^(Li,i)/Li = 
^{L2,t)/ L2. At this temperature we have, apart from 
corrections to scaling: 



A{Li,t*) ^ /Li 
A{L2,t*) \L2 



aj V 



(18) 



Note that the crossing temperature T*{Li;L2) ap- 
proaches the critical temperature for large L much faster 
than the peak of any susceptibility: t* — T*(Li; L2) — 
oc L^^^^/^ (uj is the leading correction-to-scaling 
exponent). 

From the definition [Eq. p^ ] of the correlation length, 
one sees that, respectively, S^/L ^ 0{L'^'^) in the "or- 
dered" (low T, low H) phase and ^/L - 0{1/L) in the 
"disordered" phase. The constant c is 1/2 in the case 
when the ordered phase has a Z2 global symmetry, for in 
finite lattices the disconnected susceptibility would van- 
ish. Near a second-order transition, ^/L does not depend 
on L, so there is a region in which ^/L for different lattice 
sizes cross each other. The method then consist in find- 
ing the value T*{Li;L2) at which this crossing happens 
and extracting the exponent a/v by means of Eq. psp . 

We apply the methods to several observables to extract 
exponents, in particular 



Ac 1 Ac 
-''^dis ' Xdis 
Cv 

|Afs| 



a = 7 = ^(2-7?), (19) 

a = 7 = t.(2 - 77), (20) 

a = a, (21) 

a = -P, (22) 

a = l + v. (23) 



Finite-size effects in first-order phase transitions on 
pure systems are qualitatively similar to their second- 
order counterpart, provided that one considers effective 
critical exponentsj ^^'^^i^^ With the assumption that the 
lattice size is much larger than the correlation length, 
simple scaling relations hold for the size of the broadened 
transition region, the height of the peak of the specific 
heat and the extremal point of the binder cumulant for 
the energy density: denoting with subscripts -I- and — 
values of observables of the two competing phases at a 
first-order transition in the infinite volume limit (one of 
the phases can be degenerate) and being Q = — E- 
the latent heat, one has the following in a finite system>^ 



T*{L)^T, 
Cu{T*) 
l-gf{T*) 



a{Q)L-\ (24) 

cl(a+,a-) + c2(g)L^ (25) 

9i{E+,E_) (26) 



where, in particular, a(Q), C2{Q) and _gi(-B+,i?_) vanish 
if the latent heat is zero, i.e., E+ = are the 

specific heats of the ± phases. Finally, the susceptibility 
also diverges with the volume of the system. 

However, in the presence of disorder the scaling law of 
T* (L) — Tc should be modified (see the appendix) 



T*(L)-T, = 5(0)L-'*/2 



(27) 



This follows, for instance, from a simple mean-field ar- 
gument)^ that yields a linear relation between the criti- 
cal temperature and the number of spins in the samples. 
Since the average spin density fiuctuates as L"''/^, we ex- 
pect this to be the width of the critical region on finite lat- 
tices. Furthermore, the specific heat and the connected 
susceptibility may diverge only as fast as L'^''^. See the 
appendix for a detailed discussion of these bounds. 

Hence, assuming that the observables diverge as much 
as possible, we can write the following "apparent" critical 
exponents for a disordered first-order transition: 
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a 

V 

1 

V 



d 

2 ' 
d 

2 ' 
d 

2 ■ 



(28) 



(29) 



(30) 



Notice that we follow Ogielski and Husei in defining ?/. 



From the last equation and using rj = 2 — ^ jv ~ 2 — d/2, 
we find in d = 3 that ij = 0.5 and v ~ 2/3. 

If we assume that the averaged probability distribu- 
tion of the energy {P{E)) is composed (in the thermo- 
dynamic limit and at the critical point) bof the sum of 
Dirac deltas, we should obtain a divergence L"^ for the 
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FIG. 1: The connected susceptibility computed with \Ma\ as 
function of T for lattice sizes L = 8, 12, 16, 20 and 24. Lines 
are interpolating splines as a guide to the eye. 

normalized variance of this averaged probability, obtain- 
ing (e.g., for the energy) 

L"^ (W) - M') = Q^L^ ■ (31) 

In particular, we should recover Eq. (j26[) for gf which 
is computed with {P{E)) [see Eq. ([6])]. Please note that 
the width of {P{E)) is not related to the specific heat, 
which is rather related to L"^ ({E^) - {E)A . 



III. SIMULATION DETAILS 

We simulate the model using the usual Metropolis algo- 
rithm with sequential spin flip schedule and the parallel 
tempering technique.— We restricted our simulation to 
the 

H^l.hT (32) 

diagonal in the (T, H) plane in order to keep away from 
the crossover to the zero field case. This should also avoid 
problems with an oblique crossing of the transition line 
and will help in the comparison with previous numerical 
simulations.— 

The critical temperature on this diagonal stays around 
Tc = 1.5, and we simulated Nt = 29 temperatures for 
every lattice size at equally spaced (T, H) values along 
this diagonal. For smaller lattices (L = 8, 12, and 16) 
the T values were in the range [1.3,2.7], while for the 
larger sizes [L ~ 20 and 24) the temperature range was 
[1.3,2.0]. For every lattice size, 1280 samples (different 
realizations of the disorder) were simulated. Statistics is 
also doubled as our program processes two real replicas 
per sample, with the same disorder, at each (T, H) value. 

The use of optimized asynchronous multispin coded 
update routines in our programs allowed us to thermal- 
ize systems on lattices with size up to L = 24. The 
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FIG. 2: Top: the susceptibility peak height as function of 
U^^" . Bottom: Peak position as function of L"^''". Solid 
lines are the fitting power-law functions. See section HV Al for 
more details. 
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TABLE I: Parameters characterizing the simulation. See text 
for details. 

program can perform Metropolis update at a 1.3 ns/spin 
rate on a conventional 64 bits Intel CPU at 3.4 GHz. 
Of course, the use of parallel tempering (PT) slows down 
the performance of the multispin code simulation, but we 
can limit the loss of performance if we let the program 
perform PT swaps every many Metropolis lattice sweeps. 
We verified that a PT swap trial every 20 Metropolis lat- 
tice sweeps also allows for hotter replicas to decorrelate 
before the exchange with colder ones, at the cost of a 
factor of 1.5 in performances. Cluster update algorithms 
did not prove convenient due to a dramatic increase in 
total computational load. In what follows, we consider 
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simulation time units such that 20 Monte Carlo (MC) 
steps are 20 Metropolis full lattice updates plus 1 PT 
step. 

Simulating 1280 samples for the largest lattice (L = 24 
and 24 X 10^ MC steps) took about four weeks and 20 
computation nodes on the Linux cluster at BIFI. By mon- 
itoring nonlocal observables (like the susceptibilities), we 
have checked that the runs are thermalizcd: we have 
reached a plateau in all the nonlocal observables we are 
measuring. In particular, for each lattice size, let fgim be 
the total time in MC steps devoted to simulate a sample, 
the time needed to achieve equilibrium always resulted 
shorter than isini/2. Indeed, we discarded measures at 
all times t < ts\in/2. Simulation parameters are summa- 
rized in Table m 

Two further thcrmalization test were provided by the 
parallel tempering statistics: (1) we have checked that 
the temperature samples perform all the road from higher 
to lower temperatures and come back; (2) the temper- 
ature samples have stayed essentially the same Monte 
Carlo time in all the temperatures simulated. 

IV. NUMERICAL RESULTS 
A. Second-order phase transition scenario 

We will use first the old fashioned peak method and 
turn later to the quotient method. In this way we will 
obtain complementary information. 

In Fig. [T] we show the staggered magnetization con- 
nected susceptibility Xc data. Clear peaks are present 
from which it is possible to extract information on expo- 
nents 7, V and 77. There is a lot of noise in the signal for 
Xc at low temperatures for large lattice sizes {L = 20 and 
L — 24). This is almost exclusively due to the discon- 
nected part of the connected susceptibility, which is dif- 
ficult to obtain because of metastability (see Sec. lIVC"]) . 

The exact peak position and height are located by 
means of cubic polynomial interpolation, and by using 
the standard second-order phase transition equations for 
the peak and the position of the maximum of the sus- 
ceptibility (xmax oc L''^'^ and Tc — Tlnax cx; L"^/"^), we 
obtain 

^ = (2 -,7) = 1.6(1) ^77 = 0.4(1), (33) 

V 

J/ =1.0(3), (34) 

rc = 1.58(8), (35) 

(data for L ~ 8 has been excluded in determining Tc and 
ly). These results are fully consistent with /3 = 0, even if 
we are using a second-order ansatz in the analysis. 

Figure [5] shows the dependence of peak heights and 
positions on the size L. These estimates are compatible 
with previous ones by Ogielski and Huseji Tc = 1.50(15), 



ly = 1.3(3), and 77 = 0.5(l)i^ Ground states calculations 
by Hartmann and Nowak^ for the DAFF (but at dilution 
p = 0.55) gave ly = 1.14(10). 
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FIG. 3: The specific heat as function of T for various lattice 
sizes. The main feature is in the decrease and flattening of 
the peaks as L increases. 

The specific heat shows no tendency to diverge at all 
near the transition region. On the contrary, the peak 
of C„ tends to slightly decrease and broaden as system 
size increases (see Fig. This is probably an artifact 
due to the large slope of the path in the (T, H) plane 
that we simulated [sec Eg. ([5^ ]. In fact, the larger the 
system size, the lower the peak height. However, note 
that the specific heat has a contribution (from the mag- 
netic energy) with an explicit linear dependence on the 
field strength (also the magnetization depends strongly 
on it). Anyhow, this supports a scenario of negative 
(maybe vanishing) a, as reported, for example, by Rieger 
and Youngii, Rieger—, Middleton and Fisher— in their 
simulation of the random field Ising model and in expcr- 
imentSfi We shall discuss further the specific heat in the 
following. At the time being, note that the peak position 
of Cy may be fitted to the usual power law and we find 

Tc^ 1.68(4), (36) 

which is compatible with the estimate given in Eq. (j35p . 
This fit provides no information on the v exponent (ly = 
2(2)). 

We can extract further information on several expo- 
nents by means of the QM. Figure [3] shows a clear cross- 
ing of the ratio ^ /L as function of T for different values of 
L. Data for i = 20 and L = 24 are quite noisy, again due 
to difficulties in measuring the disconnected part of G{k) 
[Eq. but still allow for locating a crossing tempera- 

ture with other curves. As expected on general grounds, 
the crossing temperatures stay well away from the po- 
sitions of the peak of Cy and Xc but lie fairly close to 
the T(,°° value that has been extrapolated from the peaks 
position. Indeed, in the absence of scaling corrections, 
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Li, L2 


Tcross 


V 


V 


a/v 


P/v 


8,12 


1.6(2) 


0.5(1) 


-1.0(1) 


0.091(6) 


0.07(6) 


8,16 


1.54(6) 


0.8(2) 


-0.99(3) 


0.07(2) 


0.05(2) 


8,20 


1.55(2) 


0.4(2) 


-0.97(1) 


0.070(1) 


0.049(5) 


8,24 


1.58(2) 


0.2(2) 


-0.94(1) 


0.083(1) 


0.06(1) 


12,16 


1 K /I \ 

1.5(1) 


1.1(6) 


-1.00(3) 


0.07(2) 


0.04(3) 


12,20 


1.53(4) 


0.1(4) 


-0.95(3) 


0.08(1) 


0.05(1) 


12,24 


1.57(4) 


0.0(3) 


-0.92(2) 


0.087(7) 


0.06(1) 


16,20 


1.55(5) 


-0.5(8) 


-0.90(4) 


0.09(2) 


0.056(4) 


16,24 


1.59(4) 


-0.4(5) 


-0.85(2) 


0.09(1) 


0.08(2) 



TABLE II: Exponents and crossing temperatures extracted 
with the QM applied to the intersection of the cumulant ^/i. 

there should be no system size dependency of the cross- 
ing temperature. Such dependency, if any, carries impor- 
tant information on scahng corrections Unfortunately, 
in our case, there is no clear systematic dependence of 
Tcross on the lattice size as, for example, for small sys- 
tems Tcross tends to shift to lower values as Li and L2 
increase, while it is sensibly shifted toward higher values 
when lattice size L = 20 or i = 24 is considered. This 
probably indicates that a crossover to first-order behavior 
is showing up. 

We show in Table HIl our results for values of exponents 
a, /3, "q and rj^ obtained from the QM. Unfortunately, we 
have not been able to measure dpS^ with enough precision 
to give a direct estimate for the thermal exponent v. 

Maybe the most striking result in Table HIl is the small- 
ness of the (3 exponent, indicating that the order param- 
eter could be discontinuous at the transition. A similar 
behavior has been found for the RFIM^ii^ and in exper- 
iments»i A very small (but definitely positive) value of 
/3 has been found also by Falicov et alM- by means of 
renormalization-group calculations for the binary RFIM. 
They also calculated magnetization curves as functions 
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FIG. 4: The cumulant ^/L as function of the temperature for 
various system sizes. Lines are interpolating splines and serve 
only as a guide to the eye. Note the noise in the curve for 
1/ = 24 in all the crossing region. 



of the temperature and field strength, showing abrupt 
jumps at the transition point. 

The value of rj agrees with the one found in Ref. 
0,77 = -1.0(3) and agrees with the smallness of the or- 
der parameter exponent and the estimate [Eq. ((M]) ] of 
1/ as the relation /? = (1 +77)2^/2 holds. Note also that, 
given the value of tj from Eq. (j33p . the Schwartz-Soffer— 
relations 2(77 — l)>r7>— 1 are satisfied as equalities 
within errors. We see that, at larger sizes, the value of rj 
decreases down even to negative values (showing a large 
error) but is always compatible with our previous esti- 
mate (at least at two standard deviations): r/ — 0.4(1). 
Also the disconnected susceptibility diverges as (since 
rj is very close to —1, because Xdis oc L^"**). As for the 
specific heat exponent, we know from the Harris crite- 
rion^ that a should be negative or zero in a disordered 
second-order transition framework. Here we report val- 
ues which are small but definitely positive. It is also true 
that if the specific heat had a cusplike singularity, our 
data would not allow to evaluate the asymptotic value 
for Ci,, and the estimates for a would be meaningless. 
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L=12« 




1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 

T 

FIG. 5: The energy Binder cumulant as function of T for all 
lattice sizes. 



B. First-order phase transition scenario 

The analysis presented above, based on the hypothe- 
sis that a second-order transition is taking place, looks 
inconclusive. This is especially clear from exponent 77, 
which lies so near to our prediction for a first-order phase 
transition: rj = 0.5. In addition this exponent in the 
Schwartz-Soffer inequality fixes the value of 77 to -1, and 
all our estimates of the 77 exponent are compatible with 
this value. Of course, this could be due to finite L correc- 
tions to scaling, but we think that the phase transition 
is truly first-order. 

We now proceed to show that our data are compatible 
with a weak first-order transition with a very large, but 
not diverging, correlation length at the transition point. 
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FIG. 6: Scaling of the minima of Binder cumulant Eq. p7|) 
(top) and of the minima positions (bottom) . Fitting functions 
are also showed. In the bottom we have fixed dg/2 to 3/2 in 
the fit. See section HVBI for more details. 



A good observable to test is the Binder cumulant for the 
total energy density)^ 



(37) 



which is usually easy to measure in simulations because 
of the good noise-to-signal ratio of the energy density. 
Notice that this Binder cumulant works directly on the 
averaged probability distribution of the energy. In both 
the disordered and ordered phases, well away from the 
transition temperature, the probability distribution of 
the energy {P{E)) tends to a single delta function in 
the thermodynamic limit, so that gf ^ 1. In the case 
of a second-order transition this is also true at Tc, while 
in the presence of a first-order transition, we have an 
energy distribution with more than one sharp peak, so 
the infinite volume limit of gf is nontrivial. Challa et 
al^ obtained the expression for the nontrivial limit and 
finite size correction to leading order in the framework 
of a double Gaussian approximation for the multipeaked 



l-gf{T*) 



9i{E+,E^) (38) 
g2{E+, E^,Cv+,Cy-)L^'^ , 
Et 1 
" " 2 ' 



{El- 



Et 



(39) 



where T* is the temperature at which the minimum 
(maximum of 1 — gf) appears and 172 is a complicated 
combination of the specific heats and energies of the in- 
finite volume coexistent states (-1-,—). The term gi is 
vanishing if the latent heat Q = E^ — E^ is zero. 

Our data for \ — gf effectively show broad peaks at 
temperatures T*{L) shifting toward lower temperature 
values as L increases (Fig. [S|). We also expect from 
Eq. (|27p that the critical temperature shift should scale 
as T*(L) -Tc cx L-'^''^. 

The matching of the data with this model is impressive, 
(see Fig. E]). A power-law fit against 1 - gf (T*(L)) = 
gi+g2L~'^<' gives 



51 = 

dg = 



0(2) X 10"^ 
3.0(1), 



(40) 
(41) 



where we excluded in the fit the L ~ 8 data [fitting also 
L — 8 data brings dg = 2.8(1) but gi, even if compatible 
with zero, has an unphysical negative value]. The ex- 
trapolated transition temperature is (assuming a power 
dg/2 = 1.5 and L > 8) 



Tc = 1.64(2). 



(42) 



with a reasonable x^/DOF — 0.6. We also recall that 
susceptibility data gave 1/j/ = 1.0(3) which is acceptable: 
an exponent 3/2 is within two standard deviations. 

However, the infinite volume limit of 1 — 
gE^X*{L = 00)] is very small, suggesting a zero la- 
tent heat for the transition. Actually, we will show 
below that one can estimate from gi the order of 
magnitude of the latent heat, which will turn out to be 
in agreement with metastability estimates. 



C. Metastability 

A very small latent heat may be very hard to detect due 
to large thermal and sample-to-sample (disorder) fluctu- 
ations. If this is the case, it should be possible to detect 
the latent heat by exploring the behavior of single sam- 
ples. 

Indeed, for our largest systems {L ~ 24), around 20% 
of the samples started to display metastability between 
a disordered, small Ms state and a large Ms state. This 
behavior was not detected on smaller systems. Further- 
more, for a large fraction of the samples, the metasta- 
bility on Ms was correlated with a metastability in the 
internal energy and in the magnetization density. This 
can be observed, for instance, in the Monte Carlo his- 
tory at temperature T = 1.5 {H = 2.25) shown in Fig. 
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[7] for a L = 24 sample. Note that the fluctuations for 
the internal energy were huge. However, if one bins 25 
consecutive Monte Carlo measurements (white squares in 
the central plot) metastability is very clear. 

We also learn from Fig. [7] that the probability distri- 
bution for the staggered magnetization shows three clear 
peaks, one for a disordered state and two for a quasisym- 
mctric ordered phase. The transition time is of the order 
of 10^ MC steps (and tunneling is probably sped up by 
our use of parallel tempering) . It is then clear that some 
of the samples may not have had enough time during the 
simulation (2.4 x 10^ MC steps) to perform enough tran- 
sitions between metastable states to give a correct value 
for the mean staggered magnetization, and this explains 
the noise we found in observables involving connected 
functions (susceptibilities, correlation length and specific 
heat). 

One can estimate the latent heat and the mean energy 
from Fig. [7] (recall L = 24): Q ~ 0.005 and E+ E_ = 
E ~ 1.36. We can introduce these values in the equation 
from gi [see Eq. (pS)) ]. For small latent heat (we write 
only the dominant term) it is possible to obtain 



1 



9i 



_Q_ 

AE^ 



(43) 



obtaining gi ^ 5 x 10^'*, only at two standard deviations 
of the gi value computed by extrapolating the Binder 
cumulant. 



V. CONCLUSIONS 

We have studied the three dimensional diluted antifer- 
romagnetic Ising model in a magnetic field using equilib- 
rium numerical techniques and analysis methods. 

We have found that the data can be described in the 
framework of a second phase transition, and obtained 
critical exponents are compatible with those obtained 
by Ogielski and HusCfi However, the critical exponent 
for the order parameter is very small, which points to a 
first-order transition. Note, however, that similarly small 
values of this critical exponent were found in ground- 
state investigations both for the DAFF— (at different di- 
lutions) and for the Gaussian RFIM^iii (these authors 
claimed that the phase transition was continous). Never- 
theless, by studying the Binder cumulant of the energy, 
we obtained clear indications on our largest lattices of 
a weak first-order phase transition. Furthermore, on our 
largest systems, a large number of samples show flip-flops 
between the ordered (quasidcgencratcd) and disordered 
phases both in the energy as well as in the order pa- 
rameter, which again is a strong evidence for the weakly 
first-order scenario. 

We remark that a complete theory of scaling in disor- 
dered first-order phase transitions (in line with that of 
Ref. [l^ for ordered systems) is still lacking. However, 
we have proposed a set of effective exponents and have 
shown that this scaling accounts for our data. 
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APPENDIX A: SCALING IN FIRST-ORDER 
PHASE TRANSITIONS IN PRESENCE OF 
DISORDER 

It is straightforward to use the Cauchy-Schwartz in- 
equality to obtain a bound on the p derivative of an ar- 
bitrary observable. This bound will hold both for first- 
and second-order phase transitions. Following the lines 
of reasoning of Ref. [2l[ we get 



d{A) 
dp 



< 



(Al) 



We recall that p is the dilution of the model, and = 
l/[p(l — p)]- Notice that this inequality holds for any 
temperature, dilution and lattice size. 



Assuming now that y (A^) is of the same order of mag- 
nitude of (A) (which is certainly the case for the internal 
energy), we translate Eq. (|Aip into a bound for the log- 
arithmic derivative: 



c^log(^) 
dp 



(A2) 



Now, the logarithmic derivative tells us about the width 
of the critical region on a finite system. For instance, 
at a given temperature and magnetic field, let p(E) be 
the spin dihition at a susceptibility peak and pc the ther- 
modynamic limit of any such quantity. We then expect 
p(V) — pc (X L^''/^. Notice that the notion of a critical 
region permits us to define an effective v exponent as 
p{L) — Pc (X L^^/" . Hence, 



2 



(A3) 



If the coexistence line has finite slope in the {T,p) plane, 
it is clear that the critical width in dilution is propor- 
tional to the critical width in temperature. A similar ar- 
gument holds for the derivative with respect to the mag- 
netic field. Thus, the logarithmic derivative with respect 
to temperature or magnetic field of A may diverge (at 
most) as fast as L''/^. So, we have found an upper bound 
for the divergences of the specific heat and connected 
susceptibilities (both are derivatives of the energy and 
magnetization respectively): they cannot diverge, with 
the lattice size, with an exponent greater than d/2. 
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Let us also remark that in footnote 7 of Ref. [2T| it 
is reported that v = 2/d for first-order transitions in the 
presence of disorder but without an explanation of this 
fact. 

Finally we will show that the Schwartz-Soffer— in- 
equality also holds in a first-order phase transition sce- 
nario. Schwartz and Soffcr show that 

where q is the momentum, h is the standard deviation of 
the magnetic field and Xc^° {q) and Xjfg (9) arc the Fourier 
transforms of the connected susceptibility and the discon- 
nected part of it respectively!^ In order to obtain the in- 



equality, we introduce the minimum momenta available 
which is of order 1/L, and use x(^) = Xc^^l'Zmin) and 
Xd\s{L) = Xdfs (^min). On the other hand, in a first-order 
phase transition, we define the effective exponents 77 and 
rj by means of the scaling of the susceptibilities at the 
critical point: x{L) ^ L'^^'' and Xdis(i) ^ L^^^K All 
together, we obtain the Schwartz-Soffer inequality: 



^>2-r?. (A5) 

Note that we have not used the criticality properties of 
the propagators. 
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FIG. 7: Monte Carlo histories of magnetic energy density 
(top), total energy density (middle) and staggered magneti- 
zation (bottom) for a typical sample of lattice size L — 24 
showing jumps between states at the transition temperature. 
The latent heat is clearer when binning 25 consecutive Monte 
Carlo measurements for the total energy (open symbols). 



